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A new method is presented for performing first-principles molecular-dynamics simulations of systems 
with variable occupancies. We adopt a matrix representation for the one-particle statistical operator 
F, to introduce a "projected" free energy functional G that depends on the Kohn-Sham orbitals 
only and that is invariant under their unitary transformations. The Liouville equation [T, H] = 
is always satisfied, guaranteeing a very efficient and stable variational minimization algorithm that 
can be extended to non-conventional entropic formulations or fictitious thermal distributions. 
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In recent years, the range of problems that can be stud- 
ied with quantitative accuracy using the methods of com- 
putational solid state physics has expanded dramatically. 
It is now possible to calculate many materials proper- 
ties with an accuracy that is often comparable to that 
of experiments. This degree of confidence is based on 
the fundamental quantum-mechanical treatment offered 
by density- functional theory (DFT) (lj , coupled with the 
availability of increasingly powerful computers and with 
the development of algorithms tuned towards optimal 
performance B. 

The application of these methods and techniques to 
metallic systems has nonetheless encountered several dif- 
ficulties, that have made progress slower than for the 
case of semiconductors and insulators. The discontinu- 
ous variation of the orbital occupancies across the Bril- 
louin zone (BZ) makes the occupation numbers rather 
ill-conditioned variables, and the self-consistent solution 
of the screening problem can suffer from several instabili- 
ties. The absence of a gap in the energy spectrum and the 
requirement of an exact diagonalization for the Hamilto- 
nian matrix everywhere in the BZ (in order to assign the 
occupation numbers) introduce "slow frequencies" in the 
evolution of the orbitals towards the ground state and 
preclude the straightforward extension to metals of al- 
gorithms which performed well for insulators. Smearing 
the Fermi surface with a finite electronic temperature H 
allows for an improved BZ sampling, but only partially 
alleviates the problems alluded to above. 

In this Letter, we introduce a new approach which 
solves many of these problems in a natural way, and 
which provides a general and efficient framework for ob- 
taining the ground state of a Kohn-Sham Hamiltonian 
at a finite electronic temperature. The typical context is 
the Mermin formulation for the Fermi-Dirac statistics [Q , 
but the method also applies when generalized entropic 
functionals are introduced ||, as it is often the case for 
metallic systems. Other applications include DFT stud- 
ies of insulators or semiconductors with thermally excited 



states U, and fractional quantum Hall states |7[. The 
language of ensemble-DFT [|| is used, and an orbital- 
based variational algorithm for the minimization to the 
ground state is developed and implemented. Dramatic 
improvements are obtained in the convergence of the en- 
ergies and especially of the Hellmann-Feynman forces. 

Within ensemble-DFT, the Hclmholtz free energy func- 
tional at a temperature T and for an 7V-representa- 
ble charge density n(r) in an external potential V is 
Ay[n(r)] = F T [n(r)] + J V(r)n(r)dr, where F T is the 
finite-temperature Mermin-Hohenberg-Kohn functional 
Q. The charge density n (r) that minimizes Ay is the 
ground-state charge density, and Ay[n (r)] is the free 
energy of the electronic system. A Kohn-Sham-like map- 
ping onto non-interacting electrons leads to a decompo- 
sition of the functional Ft into explicit terms (the non- 
interacting kinetic energy, the classical electrostatical en- 
ergy, and the entropic contribution) plus the unknown 
exchange-correlation functional -Br.xc, for which we take 
here the local density approximation Q. 

A key assumption is made by adopting a matrix rep- 
resentation fij, in the basis of the orbitals, for the one- 
particle effective statistical operator T; the charge density 
is correspondingly written as 



n(r) 



(1) 



Here the {ipi\ are orthonormal single-particle Kohn- 
Sham orbitals, the sum extends in principle over all the 
states, and the matrix is constrained to have trace 
equal to the number of electrons and eigenvalues bounded 
between and 1. We can then write in all generality the 
functional A to be minimized as 
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the Hartree and exchange-correlation terms, which de- 
pend on the charge density, have been grouped together. 
The entropic term is taken to be a function of the eigen- 
values of f; in the Fermi-Dirac case, it is S[{fij}) — 
trs(f), where s is /In/ + (1 — /)ln(l — /). In most ap- 
plications, the external potential V ex t is generated by an 
array of non-local ionic pseudopotentials. The free energy 
functional defined in Eq. ^ is in the form of traces of op- 
erators, and so it is covariant under a change of represen- 
tation (i.e., for a unitary transformation U of the orbitals 
{tpi})', this can be verified by letting f — > f = U f XJ> and 
— > \ipj)' = J2 m Ujm'tPrn- The covariance of the free 
energy functional A allows for the definition of a new, 
projected functional G, that depends only on the orbitals 

G[T: {^} ] := min A [T; {^} , {/y} ] . (3) 

G is invariant under any unitary transformation of the 
{ipi}'. the transformed orbitals cannot lead to a different 
value for G, by virtue of the covariance of A. 

The projected functional G represents a much better 
conditioned choice than the original free energy A when 
it is used in minimization algorithms that are based on 
the orbital propagation towards self-consistency, as is the 
case in the Car-Parrinello or conjugate gradients meth- 
ods The reasons are several, albeit related. (1) 
The functional G no longer depends on the occupan- 
cies of the orbitals or on their unitary transformations 
("rotations") in the occupied subspace. These are ill- 
conditioned, non-local degrees of freedom, with the added 
non-linear constraint of charge normalization. (2) The 
fij in this formalism have become dependent variables, 
implicitly defined by the minimization in Eq. (|5j), and 
this dependency does not enter into the calculation of 
the functional derivatives SG/Sip*, since the contribu- 
tions (dG/dfki)(dfki/dtfj*) are zero because of the mini- 
mum condition. (3) The occupancies of the orbitals and 
their rotations in the occupied subspace are now consis- 
tently considered as part of the same problem, namely 
that of finding the ground-state statistical operator, and 
not — as it is usually done — as two independent degrees of 
freedom. The evolution of the fij turn out to be largely 
decoupled from the problem of updating the orbitals in 
the subspace orthogonal to the occupied subspace. (4) 
The expensive and inefficient evolution for the orbital 
rotations is now shifted to the matrix f ; this implies that 
the associated slow frequencies in the evolution of A have 
now been compressed to zero by the minimization condi- 
tion. Subspace alignment j^] between subsequent orbital 
updates is also automatically enforced. 

This formulation naturally separates the evolution of 
the orbitals {tpi} from that of the fiy. the orbitals get 
updated in an outer loop that minimizes G to selfconsis- 
tency, and after every update an inner loop on the fij 
minimizes A at fixed orbitals {ipi}- In other words, the 



fij are projected via the inner loop onto the G surface, 
where f (or L, of which f is the matrix representation in 
the orbital basis) commutes with the non-selfconsistent 
Hamiltonian. The minimization of G is then freed from 
all constraints but the orthonormality of the {ipi}, and 
the rotations of the orbitals do not play any role. 

The commuting relationship of T and H can be made 
explicit from the minimum condition. Let us define 

hij = (^|T+v ext |^> , vp = mv^lm 

for the matrix elements of the Hamiltonian (the super- 
script [n] is a reminder that the potential v|^| c depends 
sclf-consistently on the charge density); the minimum 
condition that defines G implies 

8 A 5E^ C SS 
— — = = + — T— ii bij 

Sfji Sfji 6fji 

, f , <5£h*c Sn(r) 
= h » + J dr Mv)^f-~ T[s ^-»^ 

= ^+^ ll -TV(f)]y-M^ • (4) 

The constraint of charge conservation trf = N is taken 
into account with the introduction of the Lagrange mul- 
tiplier and the notation [s'(f)]y is used in place of 
dtrs(f) jdfji . The stationary condition in (^) is thus 

hij+vff -T[ S '(f)]ij = fj,Sij . (5) 

In] 

It follows that fij and the Hamiltonian hij + Vfe are 
diagonalized by the same unitary rotation, at fixed or- 
bitals, and thus represent "commuting" operators; the 
non-selfconsistent Liouville equation [L, H] = is satis- 
fied. The relation (||) does not imply that hij + v}^ and 
fij are diagonal, but just that there is a common trans- 
formation that diagonalizes both; the formalism per se is 
not linked to a preferred diagonal representation. 

The inner loop for the update of the occupation ma- 
trix f^ is carried out at fixed orbitals, and so it does 
not require the calculation of new matrix elements for 
the kinetic energy operator or the non-local pseudopoten- 
tial, and there are no orthogonalizations involved. The 
Fourier transforms of the {ipi} can also be eliminated by 
storing their real-space representation. We have chosen 
for the fij an iterative minimization that has a simple 
and appealing rationale: if the problem were not self- 
consistent, the solution for the equilibrium fij would be 
found by straightforwardly diagonalizing the Hamilto- 
nian matrix, calculating from its eigenvalues the ther- 
mal distribution of the occupation numbers, and rotating 
them back into the current orbital representation. Since 
the problem is self-consistent, this will not be the actual 
solution, but we use it as a search direction for a direct 
line minimization in the multidimensional space of the 
fij. The procedure is organized as follows. The matrix 
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hij, with the kinetic-energy and non-local contributions, 
is determined once for all before entering the inner loop. 
The updated charge density (let us assume that the m th 
iteration in the inner loop is taking place) is calculated 

as 



(r) 
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ij 



4rV*(r)^-(r) 



(G) 



The Hartree and exchange-correlation energy and 
potential V^™'(r) are then calculated, and the matrix 
representation is constructed. The entropy S^ 7 ™) is 
also computed, following a diagonalization of f , 
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EvHt f(m) v (m) 
1 n h I ij 
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Actually, in a traditional plane-wave approach the charge 
density (^|) is calculated more efficiently in this repre- 
sentation in which the fij are diagonal, since a tempo- 
rary rotation of the orbitals can then be performed on 
their more compact reciprocal-space representation. The 
Hamiltonian matrix is then updated using the new local 
terms, and diagonalized as 



H 



(m) = 



v (m) = z yl ' in e]'" 1 > Z 
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(m)f (m) y(m) 
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The non-self-consistent minimum for f would now be 



f ( m ) _ y( m )l f l\ m ) ,,\y( 



where fx is the (Fermi-Dirac) thermal distribution. We 
choose this as our search direction in the fij space, and 
a full line minimization is performed along the multi- 
dimensional segment f^ m+1 ^ = f( m ) +/3Af( m ), where 
Af M = f(m) _ f (m)_ Note that p p aram etrizes an un- 
constrained search, since trf^ = N at the end-points 
and thus, by linearity, at all j3. The search direction is 
determined by the eigenvalues of the non-self-consistent 
Hamiltonian, attracting the fij towards the representa- 
tion in which they commute with the Hamiltonian and 
towards the thermal equilibrium values that they would 
assume for this Hamiltonian. Since the search direction 
is determined by the eigenvalues/eigenvectors of H, and 
not from the occupation numbers, this current formalism 
can also be applied when generalized entropic function- 
als are defined, or when non-monotonic thermal distribu- 
tions are introduced [||. 

The direct minimization proceeds by calculating the 
free energy and its derivative along the search line at the 
two end-points /3 = and /3 = 1, taking into account the 
self-consistent variations in the charge density and thus 
in the Kohn-Sham Hamiltonian. The line derivative A' 
is Ei.A/^ (8 A/ 8 fa), where 
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A'(j3 = 0) is always smaller than 0, and so the iterative 
update of f takes place in a strictly variational fashion. 
We then calculate the charge density n^ m \r), the matrix 
elements V^j , and the free energy A corresponding to 
(3=1, together with the line derivative A'(f3 = 1) via 
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Since the kinetic energy and the pseudopotential con- 
tributions are exactly linear along the search direction 
(only the prefactors fij change), and the Hartree energy 
is quadratic, while the remaining exchange-correlation 
and entropic terms are very well-behaved, a cubic or a 
parabolic interpolation locates the value of (3 correspond- 
ing to the minimum with very high accuracy. More im- 
portantly, the choice of a direct minimization for f along 
a linear search implies that level-crossing instabilities are 
completely eliminated, even in the limit of zero tempera- 
ture and/or in the absence of an entropic term. In prac- 
tice, we find that one or two iterations in the inner loop 
are the optimal choice even for large systems (e.g., the 
diffusion of an adatom on a slab of 144 atoms p0[), since 
we also need self-consistency with respect to {^i}- 

G can be minimized very efficiently with a direct all- 
bands conjugate-gradients method; however, since it has 
a much broader spectrum for its eigenvalues than in an 
insulating case, it exhibits a slower convergence. This 
is essentially due to the occupancies of the higher bands 
being close to zero. To solve this problem, we have re- 
sorted to a preconditioning strategy: we choose a set 
of scaled variables in which the functional has a more 
compressed spectrum, and the search directions are cal- 
culated in this new metric. In the diagonal representa- 
tion the total energy around the minimum is a quadratic 
form J2 i n fiCiC^ n (ci^ n is the expansion coefficient of ipi 
in the n-th element of the basis set); if the steepest- 
descent directions are chosen according to scaled vari- 
ables Ci — \ffi c i,n, the preconditioned steepest-descent 
directions for the original variables become jnj 

*G 1 SG 

Sip* fi Sip* 

With some degree of overcorrection jfjj , these can also be 
used to construct conjugate directions. A generalization 
to our case is obtained by calculating the steepest-descent 
directions gi of G, passing them into the diagonal rep- 
resentation where they can be preconditioned as in (|lC|), 
and transforming them back in the initial representation. 
The steepest-descent directions are 

g* - ~ = ; £ = , (ii) 

where the primed term refers to the diagonal representa- 
tion (f = flfSij = Uf LP). The preconditioned steepest 
descents G' and G, are thus 
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FIG. 1. Convergence of the total free energy (upper panel, 
semi-logarithmic scale) and of the force on the surface atom 
(lower panel) in the 15-layer Al(llO) slab. Grey line, ABV; 
thin solid line, ensemble-DFT with 2 iterations in the inner 
loop; thick solid line, ensemble-DFT with 4 iterations in the 
inner loop. 



G; 



-H 



G ( = ^2u£G' n = -Hhfc 



(12) 



(13) 



Such preconditioned gradients G, greatly improve the 
convergence rate, given that higher bands are now up- 
dated with the same speed as lower filled bands, and 



are much cheaper to compute than the gj in Eq. (11). 
In addition to this occupancy preconditioning, a stan- 
dard kinetic-energy preconditioning should also be used 
in plane-wave calculations. One iteration on the orbitals 
consists thus of several operations: 1) each precondi- 
tioned gradient — H|^,) is calculated, conjugated with 
the previous search direction, and projected out of the 
subspace spanned by the orbitals to assure first-order 
orthonormality along the search; 2) the first derivative 
of the free energy along the multidimensional (all bands, 
all plane-waves, and all k-points) line is calculated, and a 
trial step along the search line is taken; 3) after reorthog- 
onalizing the orbitals, the new free energy provides the 
third constraint to identify the optimal, parabolic mini- 
mum along the search line. 

The complete algorithm provides a remarkably robust 
and efficient convergence. As a paradigmatic case we 
present here results for a unit cell that is 32 A long, 
and contains a 15-layer lxl Al(llO) slab. We use the 
single k-point ^(j, \, |), a fictitious Gaussian temper- 
ature p3| of 4 eV, and 64 orbitals. The large value of 
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FIG. 2. Conservation of the constant of motion in a molec- 
ular-dynamics run for the 15-layer Al(110) slab. The bottom 
curve is the electronic free energy; the top two curves add to 
that the kinetic energy of the ions (the grey line is for ABV, 
the solid thin line is for ensemble-DFT with 2 iterations in 
the inner loop). 



the temperature assures that the coarse sampling is suf- 
ficient; similar results are obtained with smaller and more 
physical temperatures. Fig. [j] monitors the convergence 
of the total free energy and of the Hellmann-Feynman 
forces as a function of the number of iterations on the 
{ipi}] we compare an optimal all-bands variational imple- 
mentation (ABV) |f| with the scheme that we have de- 
scribed (ensemble-DFT) jl4| . The improved convergence 
for the total energies, and particularly for the Hellmann- 
Feynman forces, is clearly apparent. It should be stressed 
that, at variance with the ABV case, the behavior of the 
total free energy in the line searches in both the outer and 
the inner loops is accurately parabolic; interpolated min- 
ima are usually fractions of a percent off their true val- 
ues. The improved convergence of the ionic forces leads 
in particular to a much tighter conservation of the con- 
stant of motion in molecular dynamics simulations. We 
show in Fig. || the results of a run for our Al(llO) slab. 
The timestep is 2 fs; the ions are moved after a fixed 
tolerance in the convergence of the free energy is reached 
(identical results are obtained if a fixed number of iter- 
ations is used). The systematic drift of the constant of 
motion stabilizes after ~ 0.3 ps of thermalization to —0.6 
eV/cell/ps for the ABV case, and to -0.0008 eV/cell/ps 
for ensemble-DFT. Such stability opens the way to inex- 
pensive molecular dynamics simulation of large metallic 
systems even on common workstations ficfl . 
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